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ABSTRACT 


Analysis of electromagnetic fields in inhomogeneous media is of 
practical interest in general scattering and propagation problems and in the 
study of lenses. For certain types of inhomogeneities, the fields may be 
represented in terms of two scalars (components of Hertzian potentials). 

In a general orthogonal coordinate system, these potentials satisfy second 
order differential equations. Exact solutions of these equations are known 
only for a few particular cases and in general, an approximate or numerical 
technique must be employed. The present work reviews and generalizes some 
of the main methods of attack of the problem. The results are presented in n 
form appropriate for numerical computation. 
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I. INTRODUCTION 

The problem of wave propagation in inhomogeneous media is important 
in many situations in which electromagnetic waves are caused to propagate 
through media whose index of refraction has arbitrary spatial variations. 

Typical examples arise in radar signal intercepting ionized wakes, radio 
waves entering the ionosphere, microwave diagnostics of ionized gas, inhomo- 
geneous ly-loaded waveguides as well as in the study of lenses. Electromagnetic 
fields in these inhomogeneous, and more generally anisotropic, media may be 
represented in terms of two Hertz vectors ^l n . Since the solution of Maxwell’s 
equations is uniquely determined by assigning initial values to the components 
of the electric field (E) and magnetic field (H) at t = 0 which in the meantime 
must satisfy V» e E = V* - 0, where n and e are the permeability and the 
permittivity, respectively, the degree of freedom in essence is reduced to four 
arbitrary functions [ 2 . Employing the gauge transformation, one may represent 
the field in terms of two components of either the electric (SM or magnetic (tt^) 

Hertz vectors or in terms of one component of i rr and another of H . The dif- 

e m 

ferential equations satisfied by these components are in general of higher order 
than the Becond, but under certain conditions they reduce to second order . 3 . 

If the medium parameters (ju and e) satisfy some further restrictions, the 
method of separation of variables may be employed to reduce the problem to 
a set of second order ordinary differential equations. The objective of the 
present work is to study the main available methods employed in order to solve 
such equations. The emphasis will be on the method of formulation rather than 
on different approximate or asymptotic techniques. 

Exact solutions of the equations encountered in dealing with inhomogeneous 
media are known for only a few particular profiles. In reference to this the 
important contributions of Heading and Westcott j_5j- [9 J are worth men- 
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tioning. Since only a few cases yield themselves to exact solutions, numerical 
or approximate methods are expected to play an important role. This enhances 
rather than diminishes the importance of those problems for which exact solu- 
tions can be found since exact solutions are often useful as starting points for 
approximate calculations. In addition, they may also help to establish limits of 
validity for various approximation methods. Even for cases when the exact 
solution is known, computation of such functions may not be the most economical 
way to find the solution. 

The two approaches used to solve the equations are the integral equation 
and the differential equation approach. Some of the baBic computational aspects 
related to these two approaches are considered next. 
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n. INTEGRAL EQUATION APPROACH 


The differential equation encountered in the present study has the form 
0" + P(x) 0’ + Q (x) 0 = 0 (1) 


which is valid for the range a < x < b, where the prime indicates derivative 
with respect to x; P(x) and Q(x) are functions of q and e. It may be pointed 

out that the term involving 0’ may be eliminated upon making the substitution 

x 

~ ~2 f P * r,)d71 - 

0 = 0 e . The function Q (x) may be written in one of die following 


forms 


r 




Q(x) = q{x) -T (x) = < 


q > - T(x) 
av < 


. q (x) 
^ ap J 


( 2 ) 


where q Q corresponds to the homogeneous case, q is the average value of 

Q(x) over the interval a <x <b and q (x) represents a convenient approxi- 

ap 

mation of Q(x). Obviously, of all these choices of q(x), the last is the most 
appropriate ns far as the rate of convergence of a perturbation technique is 
concerned. T(x) may be considered as a perturbation function. The solutions 
corresponding to T (x) - 0 are assumed to be in terms of known functions which 
are denoted by f (x) and g(x). By a proper choice of f and g and using Green's 
function, the solution of (1) and (2) takes the form 


I 

a 


0(x)= cf(x)+ \| V(s) K(x, s) 0(s)ds 


( 3 ) 


where c is an unknown constant to be determined from the boundary conditions. 
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\ (expansion parameter) and V(s) are related to the Wronskian of the homo- 
geneous equation (unperturbed) and the Kernel has the form 

r f(x) g(s) x > s 

K(x, s) = ^ (4) 

f(s) g(x) x < s 


Usually the Bora approximation is used to solve the resultant Fredholm Integral 
Equation by means of a Neumann series in \ [lO] - [ll] . This perturbation 
expansion has a limited radius of convergence determined by the lowest eigen- 
value of the homogeneous equation [l2[]. 

Following Drukarev [l2[] , if the function 


{ 


+ \ f(x) I V(s) g(s) 0(s) ds 


is added to (3) and 0 (s) is written in the form 

0(x) = N U(x) 

then U (x) satisfies the Volterra equation 


(3) 


(fi) 


/ 


r, 


U (x) = f(x) - \ I V(s) I f (x) g(s) - g(x) f (s) j U(s) ds 


n 


j 


and 


N - c 


f 


1-\| V (s) g(s) U(s) ds 


1-1 


'J 


(T) 


(8) 


If the Neumann's series is used to solve (7), the solution converges for all X 
and an expression for 0(x) may be obtained ns the ratio of two power series 
in \ 12 . Brysk [l2 has shown that Drukarev solution coincides with the 
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determinantal solution of the Fredholm equation. Thus, an Iterative {succes- 
sive approximation) solution of the Volterra equation generates the Fredholm 
determinants, just as an iterative solution of 0<x) generates (he Born series. 

In other words, the effort normally expended on the Born approximation suffices 
to obtain the Fredholm solution with the Drukarev transformation. 

It may be pointed out that a numerical solution of (7) is much simpler 
than (3) in that U(x) is expressed in terms of its value for s no greater than 
x, thus starting at x = a, the solution can be developed successively to larger 
x jj.3]. 
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IE. DIFFERENTIAL EQUATION APPROACH 


Equation (1) may be solved, using approximate or numerical techniques, 
subject to the appropriate boundary conditions which 0(x) must satisfy. The 
problem may be reduced to solving a set of coupled first order differential 
equations satisfied by the field components and an appropriate method, e, g. , 
finite difference, may be used to get the solution [mI. A formal solution to 
the problem may be found if the original equation is reduced to Hill's equation 
[I5j . In some situations, e. g. , the slab problem, step by step integration, 
starting from an initial value, of the field equation is possible [lG_] . In general, 
it is more appropriate to introduce auxiliary functions in order to solve (1). 
These auxiliary functions are the auxiliary impedance, admittance, reflection 
and phase shift functions. The purpose of this section is to cast the differential 
equations satisfied by these functions in an appropriate form which is convenient 
for computation. 

A: Auxiliary Impedance Function : 

Let 



then, using (1), v satisfies 

2 

v* + v + Pv + Q ~ o (in) 

The term including v may be eliminated and thus reducing the computational 
effort if we introduce 


Z = 


fL 

0 


4 P 


Then, Z satisfies the equation 

z . + z 2 = 1 p 2 + i p , _ Q 

i 4 


(ID 


(12) 
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B: Auxiliary Admittance Function : 

Following the same ideas as for the auxiliary impedance function, an 
admittance function may be introduced in the form 

0 1 

Y = Q-|, + - T (13) 

where Y is given by 

Y’- Y 2 = Q + |t’ +^T 2 (14) 

and 

T = P+^ (15) 

Due to the appearance of 0 in the denominator of Z, 7. - oo whenever 
0 " 0. Similarly, Y is infinite if 0' is equal to zero. In order to get around 
this difficulty, Garbacz [l?3 suggested the use of the equations satisfied by 7 
whenever Y tends to infinity and vice versa. A more convenient way is through 
the use of an auxiliary reflection function. 

C: Auxiliary Reflection Function 
Introducing the function 

R = ^070»“ (1G) 

where \p = , g 3 exp - J (Q - 1) dx , u = - , the function R then satis- 

fies 

R' + R 2 = g 2 (l + P + Q) -u’+ w 2 (17) 

Contrary to Z and Y, R is bounded as may be observed from (16), It maybe 
noted that all the auxiliary functions introduced satisfy similar Ricatti type 
equation whose right hand side is a known function of position. 
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D: Auxiliary Phase Shift Function: 

The phase shift concept extensively used in Quantum Mechanics L 18 , 
was originally introduced in dealing with potential scattering from cylindrical 
or spherical regions. If the solution is presented in the form 


0(x) = A<x) 


cos 6(x) J (kx) 
n 


+ sin 5 (x) Y (kx) 

n 


(18) 


where o(x) is unknown function, the asymptotic behavior of 0(x) is proportional 

1 7T 

to cos kr **(r> +— ) - - 5 , Hence 6(x) was called the phase shift. Integral and 
variational expressions 18] as well as upper and lower bounds of this function 
are known in the literature [l9] . Recently, Shafai [20], [2l] derived the 
differential equation satisfied by the phase shift due to scattering from cylindri- 
cal as well as spherical regions. The purpose of this section is to extend Shafai's 
work and to solve equation (1) using the phase-shift approach. 

Using Lagrange's method [22], a solution of (1) may be written in the 

form 

0 (x) - a (x) f (x) + |3 (x) g (x) (19) 


where 


f + j3' g = 0 


( 20 ) 


o T f f + JS' g' = T 


( 21 ) 


Equations (20) and (21) lead to the evaluation of o(x) and /3 (x). Introducing the 
amplitude function A (x) and the phase function 5 (x) given by 

a (x) = A (x) cos6(x) (22) 

£t(x) = A(x) sin6(x) (23) 

and using (20) and (21), one may show that 


5'(x) = 


1M 

W(f. g) 



6 f + sin 5 



(24) 
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A' (x) = -A (x) [cos 6 f + sin 6 gj [sin 6 f - cos 5 gj (25) 

where W(f, g) = f g* -gf' is the Wronskian. Hence, A(x) and 6(x) satisfy 
first order differential equation and notably the differential equation satisfied 
by 6(x) is independent of A(x). The monotonic characteristic of 6(x), as a 
function of the perturbation T(x), may be noted from equation (24). 

In order to solve equation (24) for 6 (x), an initial value is needed 
which may be evaluated from the boundary condition on 0(x). Assuming that 
such condition requires 

a0'+b0 = O at x = x^ (26) 


it follow, upon using (20), that 


tan Sfx^ = - 


af (Xj) + bffx^ 
ag* (x 1 ) +bg(x 1 ) 


(27) 


On the other hand, if g tends to infinity as x — » x^ while the field is assumed 
to be finite there, an initial value of zero may be taken for 6 (x). Once an 
initial value of 5(x) is known, equation (24) may be easily solved numerically 
using, for example, the Predictor-Corrector of the Fourth order Runge- 
Kutta method. 

Once 6 (x) is known, A (x) is then given by 


A (x) 


: A exp 
o 


4 


( - ~ | (cos <5f + sin6g)(sini5f“COs5g) 


dx 
W(f, g) 


(28) 


where A is the value of A(x) at x = x which is assumed to be known frmi 
o o 

the boundary conditions. Assume that we are dealing with two regions whose 
boundary is at x = x^ and their corresponding solutions are 
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0=a f(x) + £ g(x) 
o o 


in region I (x > x ) 
o 


(29) 


and 


0--A cos 6 f(x) + sin 5 g (x) in region n (x < x ) 


(30) 


where region I is assumed to bo homogeneous (unperturbed). If the boundary 
condition i 
(20), that 


condition at x = x requires the continuity of both 0 and 0', it follows, using 
o 


and 


a - A (x ) cos 6(x ) 
o o o 


6 * A (x ) sinS(x ) 
o o o 


(31) 


(32) 


Hence if a is known, usually related to the incident field, as well as 6(x ), 
o o 

A(x ) is obtained directly from (31). It may be noted that /3 may be written 
o ° 

in terms of 6 (x q ) only (0 Q = tan 6 (x q )) which leads to the significant con- 
clusion that only the solution for 6 (x) is needed if we are interested only in 
the field in region I. 
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IV. CONCLUSIONS 

The present work deals with some of the computational aspects of the 
two basic approaches used to solve a second order ordinary differential equa- 
tion. These approaches are, the integral and the differential equation tech- 
niques. In cases when the exact solution is unknown, the integral equation 
technique is more appropriate if an analytical expression via perturbation 
is attempted. Generally speaking, the mechanics of solving differential or 
integral equations by means of digital computers is such that it is preferable 
to solve non-linear initial value problem as opposed to linear problems with 
two-point conditions. The reason for this is that the initial value problem can 
be resolved by means of a simple iteration procedure, which is ideally suited 
to digital computers, whereas the two-point boundary problem requires the 
solution of a large system of equations. In addition, a test of the convergence 
of the solution for the initial value problem is much easier. Another advan- 
tage is that a solution may be obtained for intermediate dimensions of the range 
considered, whereas in two-point boundary formulation the solution must be 
carried out for each range. For some particular cases, e. g. , the slab problem, 
the reduction to an initial value problem is easy but in general an auxiliary 
function should be introduced to achieve this purpose. Since they are, in 
general, well behaved and bounded functions, the reflection coefficient and the 
phase shift functions are the most appropriate for numerical evaluation. If an 
initial value is known, it may be convenient to use the reflection coefficient 
function since it does not require computation of the homogeneous solutions 
f(x) and g(x) which are needed for the phase shift approach. In general, an 
initial value for 6 (x) is easier to evaluate since the equation for 6 (x) is 
independent of the amplitude function. 
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